Emotional Memory in Daily Life

1 Introduction

This notebook contains analyses from the MEDAL study, assessing memory bias in remitted and depressed individuals. The pre-registration for the main analyses, as well as a preprint version of the accompanying article can be found here as well. The main aim of the current work was to examine how three groups (control, remitted, depressed) differ in terms of emotional memory dynamics in a real-life setting.

#Load packages 
renv::activate()

* Project 'Z:/medal_membias' loaded. [renv 0.15.5]

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(janitor)

Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

    chisq.test, fisher.test

library(readxl)
library(data.table)

data.table 1.14.4 using 8 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attaching package: ‘data.table’

The following objects are masked from ‘package:dplyr’:

    between, first, last

library(ggplot2)
library(psych)

Attaching package: ‘psych’

The following objects are masked from ‘package:ggplot2’:

    %+%, alpha

library(devtools)

Loading required package: usethis

library(pastecs)

Attaching package: ‘pastecs’

The following objects are masked from ‘package:data.table’:

    first, last

The following objects are masked from ‘package:dplyr’:

    first, last

library(cli)
library(lmerTest) #linear mixed models 

Loading required package: lme4
Loading required package: Matrix

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

library(performance) #fits of residuals 
library(DescTools)

Registered S3 method overwritten by 'DescTools':
  method        from       
  print.palette wesanderson

Attaching package: ‘DescTools’

The following objects are masked from ‘package:psych’:

    AUC, ICC, SD

The following object is masked from ‘package:data.table’:

    %like%

library(tidyr) # Separate trigger column

Attaching package: ‘tidyr’

The following objects are masked from ‘package:Matrix’:

    expand, pack, unpack

The following object is masked from ‘package:pastecs’:

    extract

library(sjPlot) # for nice tables and interaction plots 

Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!

library(knitr) # to print the tables within notebook
library(foreach) #run parallel loops

Attaching package: ‘foreach’

The following object is masked from ‘package:DescTools’:

    %:%

library(parallel)
library(doParallel) #run parallel loops

Loading required package: iterators

library(ggpubr) # emmeans
library(mediation) #mediation analysis 

Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:dplyr’:

    select

Loading required package: mvtnorm
Loading required package: sandwich
mediation: Causal Mediation Analysis
Version: 4.5.0


Attaching package: ‘mediation’

The following object is masked from ‘package:psych’:

    mediate

library(plotly) #for interactive plots

Attaching package: ‘plotly’

The following object is masked from ‘package:MASS’:

    select

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(jtools) #for theme_apa

Attaching package: ‘jtools’

The following object is masked from ‘package:DescTools’:

    %nin%

library(JSmediation) #moderated mediation
library(wesanderson) 
source("functions.R")

2 Data

Here we load in our data and preprocess it in steps. We first need to separate the different surveys that were delivered to bring them into the same sampling times, merge in some descriptives we need as covariates, and estimate subject-centered and lagged variables for our models later.

1. Load and clean data

# Load the data
MEDAL_ESM_compleet <- read_excel("data/MEDAL_ESM_compleet.xlsx")
# Move variables around and clean names
MEDAL_ESM_compleet= MEDAL_ESM_compleet %>%
  rename(id='Movisens-ID') %>% 
  clean_names() %>%
  rename(subjectcode = 'deelnemersnr') %>% 
  relocate('subjectcode', .after = 'id')
df_medal_try = MEDAL_ESM_compleet

2. Split different EMA surveys we have

# Sleep EMA Qs
df_medal_sleep = df_medal_try %>% 
  filter(form=='Slaap') %>% 
  dplyr::select(c(1, 2, 4, 5, 7, 9, 14)) 
# Recent EMA Qs
df_medal_recent = df_medal_try %>% 
  filter(form=='Recent') %>% 
  dplyr::select(c(1, 2, 4, 5, 7, 9, contains('rec')))
# Remote EMA Qs
df_medal_remote = df_medal_try %>%
  filter(form=='Remote') %>%
  dplyr::select(c(1, 2, 4, 5, 7, 9, contains('rem')))
# Mastery EMA Qs
df_medal_mastery = df_medal_try %>%
  filter(form=='Mastery') %>% 
  dplyr::select(c(1, 2, 4, 5, 7, 9, contains('master')))
# Standard Qs
df_medal_standard = df_medal_try %>%
  filter(form=='Standaard')  %>%
  dplyr::select(c(1, 2, 4, 5, 7, 9, (-contains('mastery') & -contains('rem') & -contains('rec')))) %>% 
  dplyr::select(-c(7,13,14))

# merge sleep and standard
df_medal_merged = merge(df_medal_sleep, df_medal_standard, by=c('id', 'trigger_counter', 'subjectcode'), suffixes = c("_sleep", ""), all=T )

# merge to recent, remote and mastery
df_medal_merged= merge(df_medal_merged, df_medal_recent, by=c('id', 'trigger_counter', 'subjectcode' ), suffixes = c('', '_recent'), all = T )
df_medal_merged = merge(df_medal_merged, df_medal_remote, by=c('id', 'trigger_counter', 'subjectcode' ), suffixes = c('', '_remote'), all = T )
df_medal_merged = merge(df_medal_merged, df_medal_mastery, by=c('id', 'trigger_counter', 'subjectcode' ), suffixes = c('', '_mastery'), all = T )

# Fix start dates and times
df_medal_merged$form_start_date <- as.character(df_medal_merged$form_start_date)
df_medal_merged$form_start_date <- as.ITime(df_medal_merged$form_start_date)
df_medal_merged = df_medal_merged %>% 
  dplyr::mutate(form_start_date = ifelse(is.na(form_start_date), form_start_date_sleep, form_start_date))

3. Add variables on Education, Gender, and Age

# Load in previous data for demographics (also contains EMA)
load("data/IPAQ_Maeve_Workspace.RData")
df_MEDAL <- df_MEDAL[-c(192), ]

# Select relevant demo info
df_MEDAL = df_MEDAL %>% 
  dplyr::select(1, 3, 5:7) %>% 
  rename('subjectcode' = 'subject', 'education' = 'Education3groups', 'n_episodes'='NumberEpisodes' )

# merge demo to full EMA
df_medal_merged =  merge(df_medal_merged, df_MEDAL, by = c('subjectcode'), all = T ) %>%
  clean_names() 

df_medal_merged$education = as.factor(df_medal_merged$education)
df_medal_merged$gender = as.factor(df_medal_merged$gender)

df_medal_merged = df_medal_merged %>% 
  relocate(c('id', 'gender', 'age', 'education'), .after = 'subjectcode') %>% 
  relocate ('age', .after = 'gender') %>% 
  relocate ('education', .after = 'age') %>% 
  relocate ('id', .before = 'subjectcode')

4. Categorize the subjectcode into the three conditions (Depressed, Remitted, Control (Never-depressed))

df_medal_merged$condition = case_when(df_medal_merged$subjectcode >= 700 & df_medal_merged$subjectcode < 800 ~ "remitted",
                                      df_medal_merged$subjectcode >= 800 & df_medal_merged$subjectcode < 900 ~ "depressed", 
                                      df_medal_merged$subjectcode >= 900 ~ "control")
df_medal_merged = df_medal_merged %>% relocate ('condition', .after = 'subjectcode' )
df_medal_merged$condition <- as.factor(df_medal_merged$condition)

5. Calculate reaction time

df_medal_merged$rt_sleep = lubridate::as.difftime( df_medal_merged$form_finish_date_sleep - df_medal_merged$form_start_date_sleep)

6. Separate time from when trigger was sent

df_medal_merged = df_medal_merged %>% 
  tidyr::separate(trigger, c('Random', 'Time', 'Trigger_Time_1', 'Trigger_Time_2'))  %>% 
  dplyr::mutate(trigger_time = paste(Trigger_Time_1, Trigger_Time_2, sep=':')) %>%#separate time from previous trigger
  dplyr::mutate(trigger = paste(Random, Time, trigger_time, sep=' ')) #bring back 'trigger' column 

#Turn numbers into time format
df_medal_merged$trigger_time <- as.ITime(df_medal_merged$trigger_time)

7. Find the trigger number per day

#df_medal_merged$trigger_counter[is.na(df_medal_merged$trigger_time)] <- 1

df_medal_merged$trigger_number[is.na(df_medal_merged$trigger_time)] <- 1 
df_medal_merged$trigger_number[df_medal_merged$trigger_time >= (as.ITime("10:00:00"))& df_medal_merged$trigger_time < (as.ITime("12:00:00"))] <- 2
df_medal_merged$trigger_number[df_medal_merged$trigger_time >= (as.ITime("12:00:00"))& df_medal_merged$trigger_time < (as.ITime("14:00:00"))] <- 3
df_medal_merged$trigger_number[df_medal_merged$trigger_time >= (as.ITime("14:00:00"))& df_medal_merged$trigger_time < (as.ITime("16:00:00"))] <- 4
df_medal_merged$trigger_number[df_medal_merged$trigger_time >= (as.ITime("16:00:00"))& df_medal_merged$trigger_time < (as.ITime("18:00:00"))] <- 5
df_medal_merged$trigger_number[df_medal_merged$trigger_time >= (as.ITime("18:00:00"))& df_medal_merged$trigger_time < (as.ITime("20:00:00"))] <- 6
df_medal_merged$trigger_number[df_medal_merged$trigger_time >= (as.ITime("20:00:00"))& df_medal_merged$trigger_time < (as.ITime("23:00:01"))] <- 7

# df
df_medal_merged = df_medal_merged %>% 
  relocate ('trigger_number', .after = 'trigger_counter')%>% relocate('trigger_time', .before ='form_start_date') %>%
  relocate('trigger', .before ='trigger_time')

8. Final cleaning of the dataset

#Create dataset where double morning list is removed 
df_medal_merged = df_medal_merged %>% 
  dplyr::mutate(trigger_counter=ifelse(trigger_number==1 & lag(trigger_number)==1, lag(trigger_counter), trigger_counter)) %>% 
  dplyr::mutate(trigger_counter=ifelse(trigger_number==1 & lag(trigger_number)==1, lag(trigger_counter), trigger_counter)) %>%
  dplyr::mutate(trigger_counter=ifelse(trigger_number==1 & lag(trigger_number)==1, lag(trigger_counter), trigger_counter)) %>%
  dplyr::mutate(trigger_counter=ifelse(trigger_number==1 & lag(trigger_number)==1, lag(trigger_counter), trigger_counter)) %>% 
  dplyr::mutate(trigger_counter=ifelse(trigger_number==1 & lag(trigger_number)==1, lag(trigger_counter), trigger_counter))

# now we have the clean dataframe, remove duplicated coloumns
df_medal_clean = df_medal_merged %>%
  dplyr::group_by(subjectcode) %>% 
  distinct(trigger_counter, .keep_all = TRUE) %>% 
  dplyr::select (1:12, 17:60) %>% 
  dplyr::mutate(original_order = row_number()) %>% 
  relocate('original_order', .before = 'id')

#Set the order of the factor to control, remission, depression (for visual purposes)
df_medal_clean$condition <- factor(df_medal_clean$condition, levels=c('control', 'remitted', 'depressed' ), labels=c('control', 'remitted', 'depressed'))

# Create dummy variable 
df_medal_clean = df_medal_clean %>% 
  mutate(condition_dummy = case_when(condition == 'depressed'~ 3, condition == 'remitted'~ 2, condition == 'control' ~ 1))

df_medal_clean$condition_dummy = as.factor(df_medal_clean$condition_dummy)

#set weekday 
df_medal_clean$form_finish_date <- as.character(df_medal_clean$form_finish_date)
df_medal_clean$form_finish_date <- strptime(df_medal_clean$form_finish_date, format = '%Y-%m-%d %H:%M:%S')

df_medal_clean$finish_date <- as.Date(df_medal_clean$form_finish_date)
df_medal_clean$weekday <- as.integer(format(df_medal_clean$finish_date, '%w'))
df_medal_clean = df_medal_clean %>% relocate('weekday', .after = 'trigger_number')

#solve duplicate issue 
df_medal_duplicate = df_medal_clean %>% dplyr::group_by(subjectcode, weekday) %>% filter(duplicated(trigger_number))

  #since row 10 for subject 712 is a duplicate morning list
  df_medal_clean <- subset(df_medal_clean, !(subjectcode == 712 & original_order == 10))
  df_medal_clean = df_medal_clean %>% dplyr::group_by(subjectcode, weekday) %>% 
    mutate(temp = lead(trigger_number), 
           trigger_number = case_when(trigger_number == temp ~ (temp -1), TRUE ~ trigger_number)) %>%
    select(-temp)
  #Take the rows out that have 'NA'for the weekday since form wasn't finished and there is no data for PA, NA, or Memory
  df_medal_clean = df_medal_clean %>% filter(!is.na(weekday))
  df_medal_clean$original_order <- 1:nrow(df_medal_clean) 

# Create expanded df with all potential datapoints

df_medal_day = df_medal_clean %>% 
  dplyr::select('subjectcode', 'weekday', 'original_order') %>% 
  dplyr::mutate(weekday = as.numeric(weekday)) %>%
  dplyr::group_by(subjectcode) %>%
  dplyr:: distinct(weekday, .keep_all =T) %>% 
  dplyr::ungroup()

df_medal_day = df_medal_day %>% 
  slice(rep(1:n(), each = 7)) %>% 
  dplyr::mutate(trigger_number = rep(1:7, length.out=n())) %>% 
  dplyr::rename(order = original_order)


#merge with main dataframe
df_medal_clean =  merge(df_medal_day, df_medal_clean, by=c('subjectcode', 'trigger_number', 'weekday'), all= T) 
df_medal_clean = df_medal_clean %>% mutate(order = as.numeric(order))
df_medal_clean = arrange(df_medal_clean, order)

9. Centre, scale, and average Variables

used_vars = c("neg_mood", "rec_mem_neg_mood", "rem_mem_neg_mood", "pos_mood", "rec_mem_pos_mood", "rem_mem_pos_mood")

#Centering & scaling 
df_medal_clean = df_medal_clean %>% 
  dplyr::group_by(subjectcode) %>%
  # Center and scale
  dplyr::mutate( across(used_vars, ~ (.x/10), .names = "{.col}" ),
                 across(used_vars, ~ mean(.x, na.rm=T), .names = "{.col}_m" ),
                 across(used_vars, ~ (.x - mean(.x, na.rm=T)), .names = "{.col}_c" ), 
                 across(paste0(used_vars, "_c"), ~ DescTools::Winsorize(.x, na.rm = T), .names = "{.col}")) %>%  
  ungroup() %>%
  # Rescale to positive
  mutate(across(c(paste0(used_vars,"_c")), ~ abs(min(.x, na.rm = T)) + 1 + .x, .names = "{.col}s"))

10. Lag variables

# Remote 

df_medal_clean = df_medal_clean %>%
  dplyr::group_by(subjectcode) %>%
  mutate(rem_mem_pos_mood_lag = lead(rem_mem_pos_mood, n=7, order_by=subjectcode),
         rem_mem_neg_mood_lag = lead(rem_mem_neg_mood, n=7, order_by=subjectcode),
         pos_mood_lag_rem = lead(pos_mood, n=7, order_by=subjectcode),
         neg_mood_lag_rem = lead(neg_mood, n=7, order_by=subjectcode), 
         rem_mem_pos_mood_cs_lag = lead(rem_mem_pos_mood_cs, n=7, order_by=subjectcode),
         rem_mem_neg_mood_cs_lag = lead(rem_mem_neg_mood_cs, n=7, order_by=subjectcode),
         pos_mood_cs_lag_rem = lead(pos_mood_cs, n=7, order_by=subjectcode),
         neg_mood_cs_lag_rem = lead(neg_mood_cs, n=7, order_by=subjectcode),
         rem_mem_pos_mood_c_lag = lead(rem_mem_pos_mood_c, n=7, order_by=subjectcode),
         rem_mem_neg_mood_c_lag = lead(rem_mem_neg_mood_c, n=7, order_by=subjectcode),
         pos_mood_c_lag_rem = lead(pos_mood_c, n=7, order_by=subjectcode),
         neg_mood_c_lag_rem = lead(neg_mood_c, n=7, order_by=subjectcode)) %>% 
  filter(!is.na(original_order)) # filter out the NA rows for the lagging of recent memory 

# Recent
df_medal_clean = df_medal_clean %>% 
  dplyr::group_by(subjectcode) %>%
  mutate(rec_mem_pos_mood_lag = lead(rec_mem_pos_mood, n=1, order_by=subjectcode),
         rec_mem_neg_mood_lag = lead(rec_mem_neg_mood, n=1, order_by=subjectcode),
         pos_mood_lag_rec = lead(pos_mood, n=1, order_by=subjectcode),
         neg_mood_lag_rec = lead(neg_mood, n=1, order_by=subjectcode),
         rec_mem_pos_mood_cs_lag = lead(rec_mem_pos_mood_cs, n=1, order_by=subjectcode),
         rec_mem_neg_mood_cs_lag = lead(rec_mem_neg_mood_cs, n=1, order_by=subjectcode),
         pos_mood_cs_lag_rec = lead(pos_mood_cs, n=1, order_by=subjectcode),
         neg_mood_cs_lag_rec = lead(neg_mood_cs, n=1, order_by=subjectcode),
         rec_mem_pos_mood_c_lag = lead(rec_mem_pos_mood_c, n=1, order_by=subjectcode),
         rec_mem_neg_mood_c_lag = lead(rec_mem_neg_mood_c, n=1, order_by=subjectcode),
         pos_mood_c_lag_rec = lead(pos_mood_c, n=1, order_by=subjectcode),
         neg_mood_c_lag_rec = lead(neg_mood_c, n=1, order_by=subjectcode)) %>% 
  ungroup()
  

Descriptives

We first plot some general population descriptives, followed by some descriptives of the scales we use in the EMA weeks and compliance rates for both positive and negative mood.

Population

Demographics

#Create Descriptives Tables
summary_table = df_medal_clean %>% distinct(subjectcode, .keep_all = TRUE) %>% group_by(condition) %>% summarise(
  'N' = n(),
  'Age \n M +/- SD' = paste0(round(mean(age, na.rm =T), 1), ' +/- ' , round(sd(age, na.rm =T), 1)),
  #'NA Age' = sum(is.na(age)), 
  '% Education Low' = mean(education == 0, na.rm = T)*100,
  '% Education Middle' = mean(education == 1, na.rm = T)*100,
  '% Education High' = mean(education == 2, na.rm = T)*100,
  #'NA Education' = sum(is.na(education)), 
  '% Female' = mean(gender ==1, na.rm = T)*100,
  #'NA Gender' = sum(is.na(gender))
)

#Compliance rates
df_medal_compliance_individual = df_medal_clean %>% group_by(condition, subjectcode) %>% summarise('ComplianceRate' = sum(!is.na(trigger_number))/42*100)

`summarise()` has grouped output by 'condition'. You can override using the `.groups` argument.

length(which(df_medal_compliance_individual$ComplianceRate <70))

[1] 8

df_medal_compliance = df_medal_clean %>% group_by(condition) %>% summarise('% Mean Compliance' = sum(!is.na(trigger_number))/(n_distinct(subjectcode) *42)*100) 


# Create a table
summary_table = summary_table %>% left_join(df_medal_compliance, by = 'condition')
rempsyc::nice_table(summary_table)

condition	N	Age M +/- SD	% Education Low	% Education Middle	% Education High	% Female	% Mean Compliance
control	55	51.7 +/- 10.9	0.00	22.64	77.36	78.18	92.42
remitted	90	52.5 +/- 11.1	3.49	16.28	80.23	78.41	91.88
depressed	46	48.7 +/- 11.2	4.35	30.43	65.22	71.74	92.65

Tests

#create dataset with only unique subjectcode 
df_medal_distinct = df_medal_clean %>% distinct(subjectcode, .keep_all = TRUE)%>% select('subjectcode', 'condition', 'age', 'education', 'gender')

#Test for differences between groups 
  #Age
  df_medal_distinct %>% lm(age ~ condition, data=.) %>% summary()

Call:
lm(formula = age ~ condition, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.717  -6.511   2.489   7.489  18.283 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         51.6727     1.4876  34.736   <2e-16 ***
conditionremitted    0.8386     1.8963   0.442    0.659    
conditiondepressed  -2.9553     2.2042  -1.341    0.182    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.03 on 186 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.0192,    Adjusted R-squared:  0.008657 
F-statistic: 1.821 on 2 and 186 DF,  p-value: 0.1648

  
  #education
  chisq.test(df_medal_distinct$education, df_medal_distinct$condition)

Warning in stats::chisq.test(x, y, ...) :
  Chi-squared approximation may be incorrect

    Pearson's Chi-squared test

data:  df_medal_distinct$education and df_medal_distinct$condition
X-squared = 5.8234, df = 4, p-value = 0.2127

  
  #gender
  chisq.test(df_medal_distinct$gender, df_medal_distinct$condition)

    Pearson's Chi-squared test

data:  df_medal_distinct$gender and df_medal_distinct$condition
X-squared = 0.84533, df = 2, p-value = 0.6553

Plot

#Create a Graph  
  #age
  boxplot_age <- ggplot(df_medal_distinct, aes(x=condition, y=age, fill=condition)) +
    geom_boxplot() +
    labs(x = 'Condition', y = 'Age', tag = 'C' ) +
    ggtitle('Age per Condition') +
    scale_fill_manual(values = wesanderson::wes_palette(n=3, name='GrandBudapest2' ))  +
    theme_apa() 
  ggplotly(boxplot_age)

Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).

    
    #education
    plot_educ <- ggplot(df_medal_distinct, aes(x=condition, fill = education)) +
      geom_bar(position = 'dodge') +
      labs(x='Condition', y= 'Proportion', tag = 'A' ) +
      ggtitle('Education Level per Condition') +
      scale_fill_manual(values = wesanderson::wes_palette(n=3, name='GrandBudapest2' ), labels = c('Low' , 'Middle' , 'High'), name = 'Education')  +
      theme_apa() 
    ggplotly(plot_educ)

      
    #Gender
    plot_gender <- ggplot(df_medal_distinct, aes(x=condition, fill = gender)) +
      geom_bar(position = 'dodge') +
      labs(x = 'Condition', y = 'Proportion', tag = 'B' ) +
      ggtitle('Gender per Condition') +
      scale_fill_manual(values = wesanderson::wes_palette(n=3, name='GrandBudapest2' ), labels = c('Male', 'Female'), 'Gender')  +
      theme_apa() 
    ggplotly(plot_gender)

Positive Mood

Current Mood

describe.by(df_medal_clean$pos_mood, group=df_medal_clean$condition) %>% rbindlist() %>% mutate(vars=c("cont", "rem", "dep"))

Distribution

ggplot(data= df_medal_clean, aes(x=pos_mood)) + geom_histogram()+facet_grid(cols=vars(condition))

Centred Distribution

ggplot(data= df_medal_clean, aes(x=pos_mood_c)) + geom_histogram()+facet_grid(cols=vars(condition))

Recent Memory

describe.by(df_medal_clean$rec_mem_pos_mood, group=df_medal_clean$condition) %>% rbindlist() %>% mutate(vars=c("cont", "rem", "dep"))

Distribution

ggplot(data= df_medal_clean, aes(x=rec_mem_pos_mood)) + geom_histogram()+facet_grid(cols=vars(condition))

Centred Distribution

ggplot(data= df_medal_clean, aes(x=rec_mem_pos_mood_c)) + geom_histogram()+facet_grid(cols=vars(condition))

Remote Memory

describe.by(df_medal_clean$rem_mem_pos_mood, group=df_medal_clean$condition) %>% rbindlist() %>% mutate(vars=c("cont", "rem", "dep"))

Distribution

ggplot(data= df_medal_clean, aes(x=rem_mem_pos_mood)) + geom_histogram()+facet_grid(cols=vars(condition))

Centred Distribution

ggplot(data= df_medal_clean, aes(x=rem_mem_pos_mood_c)) + geom_histogram()+facet_grid(cols=vars(condition))

Negative Mood

Current Mood

describe.by(df_medal_clean$neg_mood, group=df_medal_clean$condition)%>% rbindlist() %>% mutate(vars=c("cont", "rem", "dep"))

Distribution

ggplot(data= df_medal_clean, aes(x=neg_mood)) + geom_histogram()+facet_grid(cols=vars(condition))

Centred Distribution

ggplot(data= df_medal_clean, aes(x=neg_mood_c)) + geom_histogram()+facet_grid(cols=vars(condition))
ggplot(data= df_medal_clean, aes(x=neg_mood_cs)) + geom_histogram()+facet_grid(cols=vars(condition))

Recent Memory

describe.by(df_medal_clean$rec_mem_neg_mood, group=df_medal_clean$condition) %>% rbindlist() %>% mutate(vars=c("cont", "rem", "dep"))

Distribution

ggplot(data= df_medal_clean, aes(x=rec_mem_neg_mood)) + geom_histogram()+facet_grid(cols=vars(condition))

Centred Distribution

ggplot(data= df_medal_clean, aes(x=rec_mem_neg_mood_c)) + geom_histogram()+facet_grid(cols=vars(condition))

Remote Memory

describe.by(df_medal_clean$rem_mem_neg_mood, group=df_medal_clean$condition) %>% rbindlist() %>% mutate(vars=c("cont", "rem", "dep"))

Distribution

ggplot(data= df_medal_clean, aes(x=rem_mem_neg_mood)) + geom_histogram()+facet_grid(cols=vars(condition))

Centred Distribution

ggplot(data= df_medal_clean, aes(x=rem_mem_neg_mood_c)) + geom_histogram()+facet_grid(cols=vars(condition))

3 Main Effects

A simple linear regression analysis was done to look at the differences in positive and negative mood ratings per condition.

3.1 Linear Mixed Models for PA and NA

#positive mood 
positive_model = lmer(pos_mood ~ condition + (1|subjectcode), data = df_medal_clean)


#negative mood
negative_model = lmer(neg_mood ~ condition + (1|subjectcode), data = df_medal_clean)



asis_output(tab_model (positive_model, negative_model, 
                      show.se = T, show.df = T, show.aic = T, transform = NULL,  
                      show.stat = T, show.std = T, 
                      title = 'Condition Differences Mood Rating', dv.labels = c('PA Scaled',  'NA Scaled') )$knitr)

Plot

#create boxplot for negative and positive affect for each group with significance levels 
combine_data = df_medal_clean %>%
  tidyr::pivot_longer(cols=c(pos_mood, neg_mood), names_to = 'mood_type', values_to = 'mood_rating') %>% 
  dplyr::mutate(mood_type = factor(mood_type, levels = c('pos_mood' , 'neg_mood')))

# Plot
lm_plot = ggplot(combine_data, aes(x=condition, y = mood_rating, fill= condition)) + 
  geom_boxplot() +
  labs(x ='Group', y ='Mood (scaled units)' ) +
  scale_x_discrete(labels = c('control' = 'Never Depressed', 'remitted' = 'Remitted', 'depressed' = 'Depressed' )) +
  facet_wrap( ~ mood_type, ncol =3, nrow = 2, labeller  = labeller(mood_type = c('pos_mood' = 'Positive Mood', 'neg_mood'  = 'Negative Mood'))) +
  geom_signif(comparisons = list(c('control', 'remitted'), c('control' , 'depressed'), c('remitted' , 'depressed' )), map_signif_level = T, y_position = c(15, 13, 11)) +
  scale_fill_manual(values = wesanderson::wes_palette(n=3, name='GrandBudapest2' ))  +
  theme_apa() + 
  guides(fill=F)
ggplotly(lm_plot)
ggsave('figure_1_main_effects.pdf', dpi = 320, width = 12, height = 8, path = "figures/")

Follow-Up

emmeans::emmeans(positive_model, pairwise ~ condition, pbkrtest.limit = 7800, lmerTest.limit=7800)
emmeans::emmeans(negative_model, pairwise ~ condition, pbkrtest.limit = 7800, lmerTest.limit=7800)

4 Hypothesis 1: Independent Models

In this section we test the first models from hypothesis 1 as stated in the pre-registraiton. These analyses look at the relationship between recent and remote emotional memory for both positive and negative affect, while also looking at group differences in these relaitonships based on depression status.

# set model families
model_families = c('gauss-link', "gauss-log", 'gauss-inv', 'Gamma-link', 'Gamma-log')
# detect cores for later
ncores = detectCores()-1

Positive Affect

First we run the positive affect models. For each of the recent and remote models, we check the best fitting model as stated in the pre-registration based on the AIC and the residuals. We then check whether mediation analyses are warranted, and run those. Within each tab below, we present the results and steps taken, as well as post-hoc analyses to examine the directionality of the effects.

Recent

# Model equation
h1_post_recent_eq = "pos_mood_cs ~ 1 + condition*rec_mem_pos_mood_cs_lag + gender + age + education + (1 + rec_mem_pos_mood_cs_lag|subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h1_pos_recent_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'), .errorhandling = "remove") %dopar% {
  fit_all_mods(h1_post_recent_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h1_pos_recent_models,  dv.labels = names(h1_pos_recent_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

summary(h1_pos_recent_models[[2]])
plot_diagnostics(h1_pos_recent_models)
car::vif(h1_pos_recent_models[[2]]) #Vifs were inspected without interaction and deemed OK

Follow-up

In the follow-up we look at the pairwise differences in the slopes, and the indiivuds

emmeans::emtrends(h1_pos_recent_models[[2]], pairwise ~ condition, var='rec_mem_pos_mood_cs_lag', lmerTest.limit = 3500, pbkrtest.limit = 3500 )

#create separate dataframes 
df_medal_clean_control = df_medal_clean %>% filter(condition == 'control')
df_medal_clean_remitted = df_medal_clean %>% filter(condition == 'remitted')
df_medal_clean_depressed = df_medal_clean %>% filter(condition == 'depressed') 


m1_pos_rec_control <-   lmer(pos_mood_cs ~ rec_mem_pos_mood_cs_lag + gender + age + education 
                     + (1  +  rec_mem_pos_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_control,
                   control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_pos_rec_remitted <-   lmer(pos_mood_cs ~ rec_mem_pos_mood_cs_lag + gender + age + education 
                     + (1 + rec_mem_pos_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_remitted,
                   control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_pos_rec_depressed <-   lmer(pos_mood_cs ~ rec_mem_pos_mood_cs_lag + gender + age + education 
                     + (1  + rec_mem_pos_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_depressed,
                   control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

asis_output(tab_model(m1_pos_rec_control, m1_pos_rec_remitted, m1_pos_rec_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

Remote

# Model equation
h1_pos_remote_eq = "pos_mood_cs ~ 1 +  condition*rem_mem_pos_mood_cs_lag + gender + age + education + (1  + rem_mem_pos_mood_cs_lag|subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h1_pos_remote_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'),  .errorhandling = "remove") %dopar% {
  fit_all_mods(h1_pos_remote_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h1_pos_remote_models, dv.labels = names(h1_pos_remote_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

summary(h1_pos_remote_models[[2]])
plot_diagnostics(h1_pos_remote_models)
car::vif(h1_pos_remote_models[[2]]) #VIFs without interactions were fine

Follow-up

emmeans::emtrends(h1_pos_remote_models[[4]], pairwise ~ condition, var='rem_mem_pos_mood_cs_lag' )

m1_pos_rem_control <- lmer(pos_mood_cs ~ 1 + rem_mem_pos_mood_cs_lag + gender + age + education + (1  + rem_mem_pos_mood_cs_lag|subjectcode), 
                     data=df_medal_clean_control,
                    control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_pos_rem_remitted <- lmer(pos_mood_cs ~ 1 + rem_mem_pos_mood_cs_lag + gender + age + education + (1  + rem_mem_pos_mood_cs_lag|subjectcode), 
                     data=df_medal_clean_remitted,
                     control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_pos_rem_depressed <- lmer(pos_mood_cs ~ 1 + rem_mem_pos_mood_cs_lag + gender + age + education + (1  + rem_mem_pos_mood_cs_lag|subjectcode), 
                       data=df_medal_clean_depressed,
                       control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

asis_output(tab_model(m1_pos_rem_control, m1_pos_rem_remitted, m1_pos_rem_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

Negative Affect

We next run the negative affect models. For each of the recent and remote models, we check the best fitting model as stated in the pre-registration based on the AIC and the residuals. We then check whether mediation analyses are warranted, and run those. Within each tab below, we present the results and steps taken, as well as post-hoc analyses to examine the directionality of the effects

Recent

# Model equation
h1_neg_recent_eq = "neg_mood_cs ~ condition*rec_mem_neg_mood_cs_lag  + gender + age + education + (1 + rec_mem_neg_mood_cs_lag | subjectcode)"
# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h1_neg_recent_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'), .errorhandling = 'remove') %dopar% {
  fit_all_mods(h1_neg_recent_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h1_neg_recent_models,  dv.labels = names(h1_neg_recent_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

summary(h1_neg_recent_models[[4]])
plot_diagnostics(h1_neg_recent_models)
car::vif(h1_neg_recent_models[[4]])

Follow-up

emmeans::emtrends(h1_neg_recent_models[[4]], pairwise ~ condition, var='rec_mem_neg_mood_cs_lag' )

emmeans::emmeans(h1_neg_recent_models[[4]], identity ~ rec_mem_neg_mood_cs_lag | condition )

m1_neg_rec_control <- glmer(neg_mood_cs ~ 1 + rec_mem_neg_mood_cs_lag + gender + age + education + (0  + rec_mem_neg_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_control,
                   family = Gamma(link=identity),
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_neg_rec_remitted <-glmer(neg_mood_cs ~ 1 + rec_mem_neg_mood_cs_lag + gender + age + education + (0  + rec_mem_neg_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_remitted,
                   family = Gamma(link=identity),
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_neg_rec_depressed <-glmer(neg_mood_cs ~ 1 + rec_mem_neg_mood_cs_lag + gender + age + education + (0  + rec_mem_neg_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_depressed,
                   family = Gamma(link=identity),
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

asis_output(tab_model(m1_neg_rec_control, m1_neg_rec_remitted, m1_neg_rec_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

Remote

# Model equation
h1_neg_remote_eq = "neg_mood_cs ~ condition*rem_mem_neg_mood_cs_lag  + gender + age + education + (0  + rem_mem_neg_mood_cs_lag|subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h1_neg_remote_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'),  .errorhandling = "remove") %dopar% {
  fit_all_mods(h1_neg_remote_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h1_neg_remote_models,   dv.labels = names(h1_neg_remote_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

summary(h1_neg_remote_models[[2]])
plot_diagnostics(h1_neg_remote_models)
car::vif(h1_neg_remote_models[[2]])

Follow-up

emmeans::emtrends(h1_neg_remote_models[[2]], pairwise ~ condition, var='rem_mem_neg_mood_cs_lag' )

m1_neg_rem_control <-  glmer(neg_mood_cs ~ 1 + rem_mem_neg_mood_cs_lag + gender + age + education + (0  + rem_mem_neg_mood_cs_lag|subjectcode), 
                   data=df_medal_clean_control,
                   family = Gamma(link=identity),
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_neg_rem_remitted <- glmer(neg_mood_cs ~ 1 + rem_mem_neg_mood_cs_lag + gender + age + education + (0  + rem_mem_neg_mood_cs_lag|subjectcode),
                   data=df_medal_clean_remitted,
                   family = Gamma(link=identity),
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m1_neg_rem_depressed <-  glmer(neg_mood_cs ~ 1 + rem_mem_neg_mood_cs_lag + gender + age + education + (0  + rem_mem_neg_mood_cs_lag|subjectcode),
                   data=df_medal_clean_depressed,
                   family = Gamma(link=identity),
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

asis_output(tab_model(m1_neg_rem_control, m1_neg_rem_remitted, m1_neg_rem_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

5 Hypothesis 2: Current Affect Moderation

For H2, we want to see whether current affective states moderate the ability to recall previous affective states. That is, do current negative affect feelings moderate the relationship between affect and recall?

#create separate dataframes for each group
df_medal_clean_control = df_medal_clean %>% filter(condition=='control')
df_medal_clean_remitted = df_medal_clean %>% filter(condition=='remitted')
df_medal_clean_depressed = df_medal_clean %>% filter(condition=='depressed')

#create a new dataframe for the remote variables, where the participants with the singular datapoints are removed 
df_medal_clean_remote = df_medal_clean %>% filter (!(subjectcode == 822 | subjectcode == 906 | subjectcode == 717))

#create separate dataframes per condition 
df_medal_clean_remote_control = df_medal_clean_remote %>% filter(condition=='control')
df_medal_clean_remote_remitted = df_medal_clean_remote %>% filter(condition=='remitted')
df_medal_clean_remote_depressed = df_medal_clean_remote %>% filter(condition=='depressed')

Positive Affect

Recent

An lmer model was fitted with gaussian identity, log & inverse as well as Gamma log, inverse, and identity. The models, however, were not a good fit and the gaussian identity model coped with singularity. Hence, a normal linear regression model was fitted.

h2_pos_recent_eq = "pos_mood_cs ~ condition*rec_mem_pos_mood_cs_lag*pos_mood_cs_lag_rec + gender + age + education + (1 + rec_mem_pos_mood_cs_lag + pos_mood_cs_lag_rec|subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h2_pos_recent_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'),  .errorhandling = "remove") %dopar% {
  fit_all_mods(h2_pos_recent_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h2_pos_recent_models,  dv.labels = names(h2_pos_recent_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

Due to singularity this was determined to be the best fit for the moddel

summary(h2_pos_recent_models[[2]])
plot_diagnostics(h2_pos_recent_models)
car::vif(h2_pos_recent_models[[2]])

Follow-up

m2_pos_rec_control <- lm(pos_mood_cs ~ 1 + rec_mem_pos_mood_cs_lag*pos_mood_cs_lag_rec + gender + age + education,
                   data=df_medal_clean_control)

m2_pos_rec_remitted <- lm(pos_mood_cs ~ 1 + rec_mem_pos_mood_cs_lag*pos_mood_cs_lag_rec + gender + age + education,
                   data=df_medal_clean_remitted)

m2_pos_rec_depressed <-  lm(pos_mood_cs ~ 1 + rec_mem_pos_mood_cs_lag*pos_mood_cs_lag_rec + gender + age + education,
                   data=df_medal_clean_depressed)

asis_output(tab_model(m2_pos_rec_control, m2_pos_rec_remitted, m2_pos_rec_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr)

Remote

An lmer model was fitted with gaussian identity, log & inverse as well as Gamma log, inverse, and identity. The models, however, were not a good fit and the gaussian identity model coped with singularity. Hence, a normal linear regression model was fitted.

h2_pos_remote_eq = "pos_mood_cs ~ condition*rem_mem_pos_mood_cs_lag*pos_mood_cs_lag_rem + gender + age + education + (1 | subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h2_pos_remote_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'),  .errorhandling = "remove") %dopar% {
  fit_all_mods(h2_pos_remote_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h2_pos_remote_models,  dv.labels = names(h2_pos_remote_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

Due to singularity a regular linear model was fitted.

summary(h2_pos_remote_models[[2]])
plot_diagnostics(h2_pos_remote_models)
car::vif(h2_pos_remote_models[[2]])

Follow-up

m2_pos_rem_control <- lm(pos_mood_cs ~ rem_mem_pos_mood_cs_lag*pos_mood_cs_lag_rem  + gender + age + education, 
                   data=df_medal_clean_remote_control)

m2_pos_rem_remitted <- lm(pos_mood_cs ~ rem_mem_pos_mood_cs_lag*pos_mood_cs_lag_rem  + gender + age + education,  
                   data=df_medal_clean_remote_remitted)

m2_pos_rem_depressed <-  lm(pos_mood_cs ~ rem_mem_pos_mood_cs_lag*pos_mood_cs_lag_rem  + gender + age + education, 
                   data=df_medal_clean_remote_depressed)

asis_output(tab_model(m2_pos_rem_control, m2_pos_rem_remitted, m2_pos_rem_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

Negative Affect

Recent

h2_neg_recent_eq = "neg_mood_cs ~ condition*rec_mem_neg_mood_cs_lag*neg_mood_cs_lag_rec + gender + age + education + (1 + rec_mem_neg_mood_cs_lag + neg_mood_cs_lag_rec|subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h2_neg_recent_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'),  .errorhandling = "remove") %dopar% {
  fit_all_mods(h2_neg_recent_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h2_neg_recent_models, dv.labels = names(h2_neg_recent_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

Check Final Model

summary(h2_neg_recent_models[[4]])
plot_diagnostics(h2_neg_recent_models) 
car::vif(h2_neg_recent_models[[4]])

Follow-up

m2_neg_rec_control <- glmer(neg_mood_cs ~ rec_mem_neg_mood_cs_lag*neg_mood_cs_lag_rec + 
                              gender + age + education + 
                              (1  + rec_mem_neg_mood_cs_lag | subjectcode) + (1 +  neg_mood_cs_lag_rec|subjectcode), 
                   data=df_medal_clean_control,
                   family = Gamma(link=identity), 
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m2_neg_rec_remitted <- glmer(neg_mood_cs ~ rec_mem_neg_mood_cs_lag*neg_mood_cs_lag_rec + gender + age + education + 
                              (1  + rec_mem_neg_mood_cs_lag | subjectcode) + (1 +  neg_mood_cs_lag_rec|subjectcode), 
                   data=df_medal_clean_remitted,
                   family = Gamma(link=identity), 
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m2_neg_rec_depressed <-  glmer(neg_mood_cs ~ rec_mem_neg_mood_cs_lag*neg_mood_cs_lag_rec + gender + age + education + 
                              (1  + rec_mem_neg_mood_cs_lag | subjectcode) + (1 +  neg_mood_cs_lag_rec|subjectcode),  
                   data=df_medal_clean_depressed,
                   family = Gamma(link=identity), 
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

asis_output(tab_model(m2_neg_rec_control, m2_neg_rec_remitted, m2_neg_rec_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

Remote

h2_neg_remote_eq = "neg_mood_cs ~ condition*rem_mem_neg_mood_cs_lag*neg_mood_cs_lag_rem + gender + age + education + (1 + rem_mem_neg_mood_cs_lag + neg_mood_cs_lag_rem|subjectcode)"

# set clusters
cl = makeCluster(ncores)
registerDoParallel(cl)

# run parallel
h2_neg_remote_models = foreach(family = model_families, .combine = 'c', .packages = c('lmerTest'),  .errorhandling = "remove") %dopar% {
  fit_all_mods(h2_neg_remote_eq, df_medal_clean, family) }
stopCluster(cl)

# print the models into a single table
asis_output( tab_model(h2_neg_remote_models, dv.labels = names(h2_neg_remote_models), 
                       show.icc = T, show.aic = T, transform = NULL, show.se = T, digits = 3)$knitr)

Check Final Model

summary(h2_neg_remote_models[[4]])
plot_diagnostics(h2_neg_remote_models)
car::vif(h2_neg_remote_models[[4]])

Follow-up

m2_neg_rem_control <- glmer(neg_mood_cs ~ rem_mem_neg_mood_cs_lag*neg_mood_cs_lag_rem  + gender + age + education
                            + (1  +  rem_mem_neg_mood_cs_lag|subjectcode) + (1 + neg_mood_cs_lag_rem|subjectcode), 
                   data=df_medal_clean_remote_control,
                   family = Gamma(link=identity), 
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m2_neg_rem_remitted <- glmer(neg_mood_cs ~ rem_mem_neg_mood_cs_lag*neg_mood_cs_lag_rem  + gender + age + education
                            + (1  +  rem_mem_neg_mood_cs_lag|subjectcode) + (1 + neg_mood_cs_lag_rem|subjectcode),
                   data=df_medal_clean_remote_remitted,
                   family = Gamma(link=identity), 
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

m2_neg_rem_depressed <-  glmer(neg_mood_cs ~ rem_mem_neg_mood_cs_lag*neg_mood_cs_lag_rem  + gender + age + education
                            + (1  +  rem_mem_neg_mood_cs_lag|subjectcode) + (1 + neg_mood_cs_lag_rem|subjectcode), 
                   data=df_medal_clean_remote_depressed,
                   family = Gamma(link=identity), 
                   control = glmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)))

asis_output(tab_model(m2_neg_rem_control, m2_neg_rem_remitted, m2_neg_rem_depressed, show.se = T, show.aic = T, transform = NULL, dv.labels = c('control', 'remitted', 'depressed'))$knitr) 

6 Exploratory Analyses

6.1 Mediation Models

Following the tests we ran for hypothesis 1, we would also like to see whether or not mediation effects exist.

RecPA

#create a new dataframe for the remote variables, where the participants with the singular datapoints are removed 
df_medal_moderation =  df_medal_clean %>% 
  filter (!(subjectcode == 822 | subjectcode == 906 | subjectcode == 717)) %>% 
  filter(!(is.na(pos_mood) | is.na(rec_mem_pos_mood_lag))) 

# med model
med_fit_h1_pos_rec <- lme4::lmer(rec_mem_pos_mood_lag ~ condition_dummy +
                                   age + gender + education + 
                                   (1 |subjectcode),
                                 data = df_medal_moderation)
# full model
out_fit_h1_pos_rec <-lme4::lmer(pos_mood ~ condition_dummy + rec_mem_pos_mood_lag +
                                  age + gender +education +
                                  (1 + rec_mem_pos_mood_lag |subjectcode), 
                                data = df_medal_moderation)

#mediation analysis for control and depressed
med_pos_rec_dep <- mediation::mediate(med_fit_h1_pos_rec, out_fit_h1_pos_rec,
                                      treat = 'condition_dummy' , control.value = 1, treat.value = 3 ,  mediator = 'rec_mem_pos_mood_lag', 
                                      sims=10000)
summary(med_pos_rec_dep)

RemNA

df_medal_moderation = df_medal_clean_remote %>% filter(!(is.na(neg_mood) | is.na(rec_mem_neg_mood_lag))) 

# med model
med_fit_h1_neg_rec <- lme4::lmer(rec_mem_neg_mood_lag ~ condition_dummy +
                                   age + gender + education +  
                                   (1 |subjectcode), 
                                 data = df_medal_moderation)
# full model
out_fit_h1_neg_rec <- lme4::lmer(neg_mood ~ condition_dummy + rec_mem_neg_mood_lag +
                                  age + gender +education + 
                                  (1 + rec_mem_neg_mood_lag|subjectcode), 
                                data = df_medal_moderation)


#mediation analysis for control and remitted
summary(mediation::mediate(med_fit_h1_neg_rec, out_fit_h1_neg_rec, 
                           treat = 'condition_dummy' , mediator = 'rec_mem_neg_mood_lag', control.value = 1, treat.value = 2 ,   
                           sims=10000))

#mediation analysis for control and depressed
summary(mediation::mediate(med_fit_h1_neg_rec, out_fit_h1_neg_rec, 
                           treat = 'condition_dummy' ,  mediator = 'rec_mem_neg_mood_lag', control.value = 1, treat.value = 3 ,
                           sims=10000))

6.2 Moderated Mediation

RemPA-CurrPA

df_medal_moderation = df_medal_clean_remote %>% 
  filter(!(is.na(pos_mood) | is.na(rem_mem_pos_mood_lag))) 

#mediation first
med_fit_h2_pos_rem <- lme4::lmer(rem_mem_pos_mood_lag ~ condition_dummy + 
                                   age + gender + education + 
                                   (1 |subjectcode), 
                                 data =df_medal_moderation)
# full model
out_fit_h2_pos_rem <-lme4::lmer(pos_mood ~ condition_dummy + rem_mem_pos_mood_lag 
                                + age + gender +education +
                                  (1 + rem_mem_pos_mood_lag |subjectcode),
                                data=df_medal_moderation)


#mediation analysis for control and remitted
med_pos_rem_rem <- mediation::mediate(med_fit_h2_pos_rem, out_fit_h2_pos_rem, 
                                      treat = 'condition_dummy' , control.value = 1, treat.value = 2 ,  mediator = 'rem_mem_pos_mood_lag')
summary(med_pos_rem_rem)

#create new dataframes with only control remitted and only control depressed
df_medal_con_rem = df_medal_clean %>% filter(condition=='control' | condition == 'remitted')
df_medal_con_dep = df_medal_clean %>% filter(condition=='control' | condition == 'depressed')

#remitted group
df_medal_con_rem$condition_dummy = build_contrast(df_medal_con_rem$condition, 'control', 'remitted')
med_h2_pos_rem_rem <- mdt_moderated(data = df_medal_con_rem, IV =condition_dummy, DV = pos_mood, M = rem_mem_pos_mood_lag, Mod = pos_mood_lag_rem) 
rbindlist(med_h2_pos_rem_rem$paths, idcol = T)

RemNA - CurrNA

df_medal_moderation = df_medal_clean_remote %>% filter(!(is.na(neg_mood) | is.na(rem_mem_neg_mood_lag))) 

#mediation first
med_fit_h2_neg_rem <- lme4::lmer(rem_mem_neg_mood_lag ~ condition_dummy +
                                   age + gender + education +
                                   (1 |subjectcode),
                                 control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)),
                                 data =df_medal_moderation)

out_fit_h2_neg_rem <-lme4::lmer(neg_mood ~ condition_dummy + rem_mem_neg_mood_lag +
                                  age + gender + education +
                                  (1 + rem_mem_neg_mood_lag |subjectcode), 
                                control = lmerControl(calc.derivs = F, optimizer='bobyqa', optCtrl=list(maxfun=1e6)),
                                data=df_medal_moderation)

#mediation analysis for control and depressed
med_neg_rem_dep <- mediation::mediate(med_fit_h2_neg_rem, out_fit_h2_neg_rem,  
                                      mediator = 'rem_mem_neg_mood_lag',  treat = 'condition_dummy' , 
                                      control.value = 1, treat.value = 3 , sims=10000) 
summary(med_neg_rem_dep)

#mediation analysis for control and remitted
med_neg_rem_rem <- mediation::mediate(med_fit_h2_neg_rem, out_fit_h2_neg_rem, treat = 'condition_dummy' , control.value = 1, treat.value = 2 ,  mediator = 'rem_mem_neg_mood_lag', sims=10000) 
summary(med_neg_rem_rem)

#create new dataframes with only control remitted and only control depressed
df_medal_con_rem = df_medal_clean %>% filter(condition=='control' | condition == 'remitted')
df_medal_con_dep = df_medal_clean %>% filter(condition=='control' | condition == 'depressed')

#remitted group
df_medal_con_rem$condition_dummy = build_contrast(df_medal_con_rem$condition, 'control', 'remitted')
med_h2_neg_rem_rem <- mdt_moderated(data = df_medal_con_rem, IV = condition_dummy, DV = neg_mood, M = rem_mem_neg_mood_lag, Mod = neg_mood_lag_rem) 
rbindlist(med_h2_neg_rem_rem$paths, idcol = T)

#depressed group
df_medal_con_dep$condition_dummy = build_contrast(df_medal_con_dep$condition, 'control', 'depressed')
med_h2_neg_rem_dep <- mdt_moderated(data = df_medal_con_dep, IV = condition_dummy, DV = neg_mood, M = rem_mem_neg_mood_lag, Mod = neg_mood_lag_rem) 
rbindlist(med_h2_neg_rem_dep$paths,idcol = T)
  

6.3 Context

asis_output(tab_model( (lmer(pos_mood_cs ~ rec_mem_pos_mood_cs_lag*context_location + (1  | subjectcode), data = df_medal_clean)), 
           (lmer(pos_mood_cs ~ rem_mem_pos_mood_cs_lag*context_location + (1  | subjectcode), data = df_medal_clean)),
            (lmer(neg_mood_cs ~ rem_mem_neg_mood_cs_lag*context_location + (1  | subjectcode), data = df_medal_clean)),
             (lmer(neg_mood_cs ~ rem_mem_neg_mood_cs_lag*context_location + (1  | subjectcode), data = df_medal_clean)))$knitr)

6.4 SRET

In addition to the EMA, participants also completed the SRET. This is a lab based memory bias measure. We want to check whether our memory bias measures from real-life are linked to ones from the lab in this supplementary analysis.

library('haven')
sret_file = file.path("Z:/", "inbox", "transfer-2023-12-07-02-15-pm", 'MEDAL_pre and post quest and remote recall_workfile!!_for paper only var.sav')
sret_data = read_sav(sret_file) 
sret_data = sret_data %>% select(c(1,4)) %>% distinct()

df_medal_rec_cor1 = df_medal_clean %>%
  dplyr::select(subjectcode, pos_mood, rec_mem_pos_mood_lag) %>%
  na.omit() %>% 
  dplyr::group_by(subjectcode) %>%
  dplyr::mutate(corr =  tryCatch(cor(rec_mem_pos_mood_lag, pos_mood), warning = function(e) NA),
                model_type = "PA_REC") %>% 
  dplyr::select(subjectcode, corr, model_type) %>% 
  dplyr::distinct()  %>% ungroup()


df_medal_rec_cor2 = df_medal_clean %>%
  dplyr::select(subjectcode, neg_mood, rec_mem_neg_mood_lag) %>%
  na.omit() %>% 
  dplyr::group_by(subjectcode) %>%
  dplyr::mutate(corr =  tryCatch(cor(rec_mem_neg_mood_lag, neg_mood), warning = function(e) NA),
                model_type = "NA_REC") %>% 
  dplyr::select(subjectcode, corr, model_type) %>% 
  dplyr::distinct() %>% ungroup()


df_medal_rem_cor1 = df_medal_clean %>%
  dplyr::select(subjectcode, pos_mood, rem_mem_pos_mood_lag) %>%
  na.omit() %>% 
  dplyr::group_by(subjectcode) %>%
  dplyr::mutate(corr = tryCatch(cor(rem_mem_pos_mood_lag, pos_mood), warning = function(e) NA),
                model_type = "PA_REM") %>% 
  dplyr::select(subjectcode, corr, model_type) %>% 
  dplyr::distinct()

df_medal_rem_cor2 = df_medal_clean %>%
  dplyr::select(subjectcode, neg_mood, rem_mem_neg_mood_lag) %>%
  na.omit() %>% 
  dplyr::group_by(subjectcode) %>%
  dplyr::mutate(corr = tryCatch(cor(rem_mem_neg_mood_lag, neg_mood), warning = function(e) NA),
                model_type = "NA_REM") %>% 
  dplyr::select(subjectcode, corr, model_type) %>% 
  dplyr::distinct()


df_cors = rbindlist(list(df_medal_rec_cor1, df_medal_rec_cor2, df_medal_rem_cor1, df_medal_rem_cor2))

df_ref = merge(sret_data, df_cors, by.x = 'subject', by.y = 'subjectcode')
df_ref = df_ref %>% rename(SRET = SRET_BiasGotlib_neg) %>% dplyr::mutate(group = ifelse(subject<800, "remitted", 
                                                                    ifelse(subject>=900, "control", 'depressed')) )


for (mod_type in c("PA_REC", "PA_REM", "NA_REC", "NA_REM")){
  print(mod_type)
  print(summary(lm(SRET ~ corr*group, data = df_ref %>% filter(model_type == mod_type))))
}

6.5 Relative Bias

Here we want to check for the relative NA/PA bias. So we extact the random effects from the models to get the slopes for each subject, and then compare the NA/PA slopes. This tells us if theres a bias between negative and positive memory.

# Get individual slopes
pos_rec_ref = as.data.frame(ranef(h1_pos_recent_models[[2]])) %>% filter(term!='(Intercept)') %>% select(c(3:5))
pos_rem_ref = as.data.frame(ranef(h1_pos_remote_models[[2]])) %>% filter(term!='(Intercept)') %>% select(c(3:5))
neg_rec_ref = as.data.frame(ranef(h1_neg_recent_models[[4]])) %>% filter(term!='(Intercept)') %>% select(c(3:5))
neg_rem_ref = as.data.frame(ranef(h1_neg_remote_models[[4]])) %>% filter(term!='(Intercept)') %>% select(c(3:5))
pos_rec_ref$mod = "PA_REC"
pos_rem_ref$mod = "PA_REM"
neg_rec_ref$mod = "NA_REC"
neg_rem_ref$mod = "NA_REM"

# Put together into a new data frame
df_ref = rbindlist(list(pos_rec_ref, pos_rem_ref, neg_rec_ref, neg_rem_ref))
df_ref = df_ref %>% 
  dplyr::mutate(sub = as.numeric(as.character(grp))) %>%
  dplyr::mutate(group = ifelse(sub <= 800, "remitted", 
                               ifelse(sub>=900, "control", 'depressed')) )

df_ref_rec = df_ref %>% filter(mod == 'PA_REC' | mod=='NA_REC')
df_ref_rem = df_ref %>% filter(mod == 'PA_REM' | mod=='NA_REM')
summary(lmer(condsd ~ mod*group + (1|sub), data = df_ref_rec))
summary(lmer(condsd ~ mod*group + (1|sub), data = df_ref_rem))

7 Plots

First we set our theme and fix the data for the figures we want

# Theme
ggtheme = theme(text = element_text(size = 8), 
                panel.background = element_rect(fill = "transparent"), 
                panel.grid.major = element_line(color = "grey90"),
                panel.grid.minor = element_line(color = "grey100"),
                panel.border = element_rect(color = "grey80", fill = "transparent"))

# Estimate the SD
df_medal_clean2 = df_medal_clean %>% 
  mutate(pos_mood_sd = ifelse( (pos_mood_cs_lag_rec > (mean(pos_mood_cs_lag_rec, na.rm=T)+sd(pos_mood_cs_lag_rec, na.rm=T))), "+1SD",
                               ifelse( (pos_mood_cs_lag_rec < (mean(pos_mood_cs_lag_rec, na.rm=T)-sd(pos_mood_cs_lag_rec, na.rm =T))), "-1SD", "Mean"))) %>% 
    mutate(neg_mood_sd = ifelse( (neg_mood_cs_lag_rec > (mean(neg_mood_cs_lag_rec, na.rm=T)+sd(neg_mood_cs_lag_rec, na.rm=T))), "+1SD",
                               ifelse( (neg_mood_cs_lag_rec < (mean(neg_mood_cs_lag_rec, na.rm=T)-sd(neg_mood_cs_lag_rec, na.rm =T))), "-1SD", "Mean")))

# Fix factor levels
df_medal_clean2$Condition =  factor(df_medal_clean2$condition, levels=c('control', 'remitted', 'depressed' ), labels=c('Never-Depressed', 'Remitted', 'Depressed'))
df_medal_clean2$pos_mood_sd =  factor(df_medal_clean2$pos_mood_sd, levels=c('-1SD', 'Mean', '+1SD' ))
df_medal_clean2$neg_mood_sd =  factor(df_medal_clean2$neg_mood_sd, levels=c('-1SD', 'Mean', '+1SD' ))

7.0.1 PA-CUR-REC

# Plot
plot_pa_rec = ggplot(df_medal_clean2, aes(x = rec_mem_pos_mood_c_lag, y=pos_mood_c, color = pos_mood_sd, fill = pos_mood_sd)) +
  geom_smooth(method = 'lm', alpha = 0.15) +
  labs( x = bquote(paste(REC[PM], " (subject-centered a.u.)")), y = "PM (subject-centered a.u.)", color = bquote(CUR[PM]), fill = bquote(CUR[PM]) ) +
  ggtheme + 
  facet_grid(.~Condition) 
plot_pa_rec

7.0.2 PA-CUR-REC

# Plot
plot_na_rec = 
  ggplot(df_medal_clean2, aes(x = rec_mem_neg_mood_c_lag, y=neg_mood_c, color = neg_mood_sd, fill = neg_mood_sd)) +
  geom_smooth(method = 'lm', alpha = 0.15) +
  labs( x = bquote(paste(REC['NM'], " (subject-centered a.u.)")), y = "NM (subject-centered a.u.)", color = bquote(CUR["NM"]), fill = bquote(CUR["NM"]) )+
  ggtheme + 
  facet_grid(.~Condition) 
plot_na_rec

7.0.3 Combined Plot for pub

ggarrange(NA, plot_pa_rec, NA,  plot_na_rec,
          widths=c(0.05, 1, 0.05, 1), ncol = 2, nrow = 2,
          labels=c("A", NA, "B"))
ggsave("figures/figure_2_threewayInteraction.pdf", device = 'pdf', dpi = 320, height = 6)